Clustering Mounds 200m

Finding mound clusters

Traditionally, two mounds within 100-200m are considered and labelled by archaeologists as a necropolis. Can we determine such clusters computationally? This script shows how to find clusters using mound buffers of 100m and combining this output with Voronoi polygons. Some weakness is in linearly distributed necropoleis, which can escape attention, so clusters need to be composed from bottom up (lowest number of members to the highest number of members)

TO DO:

• eliminate duplicates (see line 77). Best suggestion: run proximity spatial query on mounds that are within 15-20m of one another as potential duplicates.
• trickle spatial duplicate elimination into the initial analyses!!!

Load 2009-2022 features data

We have a total of 1090 mounds, extinct and uncertain mound features in Yambol.

# Yambol mounds
features <- readRDS("../output_data/mnds_aggr_27Dec.rds") # dataset is most recent from 27/12/2022
Y_region <- st_read("../data/YamRegion.shp")

## Reading layer `YamRegion' from data source 
##   `C:\Users\adela\Documents\RStudio\YambolMoundAnalysis\data\YamRegion.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 431400.8 ymin: 4641684 xmax: 504849.3 ymax: 4727299
## Projected CRS: WGS 84 / UTM zone 35N

head(features)

## Simple feature collection with 6 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 441126.7 ymin: 4706447 xmax: 454596 ymax: 4710087
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 6 × 19
##    TRAP Source              Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##   <dbl> <chr>               <chr> <chr>   <chr>  <chr>   <chr>   <chr>   <chr>  
## 1  9305 <NA>                Buri… Scrub   Scrub  25      1.8     3       <NA>   
## 2  9411 Survey              Buri… Scrub   Scrub  15      1.5     2       Nodata 
## 3  9416 Legacy verification Buri… Annual… Annua… 20      2       3       Agricu…
## 4  9417 Legacy verification Buri… Pasture Pastu… 15      2       4       Looting
## 5  9409 Legacy verification Buri… Annual… Pastu… 30      2       2       Agricu…
## 6  9314 Legacy verification Buri… Pasture Pastu… 30      2       5       Looting
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## #   elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## #   TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## #   ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition, ⁵​PrincipalSourceOfImpact

features %>% group_by(Type) %>% tally()

## Simple feature collection with 11 features and 2 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 420515 ymin: 4647347 xmax: 499612.3 ymax: 4722257
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 11 × 3
##    Type                      n                                          geometry
##    <chr>                 <int>                                    <GEOMETRY [m]>
##  1 Burial Mound           1024 MULTIPOINT ((420515 4693580), (420707.1 4690851)…
##  2 Burial Mound?            20 MULTIPOINT ((462092.8 4663311), (465604.9 465350…
##  3 Extinct Burial Mound    188 MULTIPOINT ((420832.5 4689836), (421875.6 469359…
##  4 Extinct Burial Mound?     4 MULTIPOINT ((487113.7 4684048), (488727 4683746)…
##  5 Other                    75 MULTIPOINT ((426968.4 4693384), (462171.2 465546…
##  6 Surface Scatter          64 MULTIPOINT ((421045.4 4690265), (424207.4 469496…
##  7 Tell                      3 MULTIPOINT ((465676.3 4675875), (467848.2 466240…
##  8 Tell?                     1                          POINT (483400.7 4697154)
##  9 Uncertain Feature        36 MULTIPOINT ((437475.3 4690926), (445710.4 469931…
## 10 Uncertain Feature?        1                          POINT (480327.6 4699085)
## 11 Uncertain Mound          68 MULTIPOINT ((422391.7 4698832), (429156.7 469547…

Clean out uncertainties

features <- features %>%
  #st_drop_geometry() %>%  
  dplyr::mutate(Condition = str_extract(Condition, "\\d")) %>%
  dplyr::mutate(Condition = case_when(Condition == 0 ~ "NA",
                               Condition == 6 ~ "5",
                               Condition != 0 ~ Condition)) %>% 
  dplyr::mutate(Condition = as.factor(Condition)) %>% 
  dplyr::mutate(Type= gsub("\\?","",Type)) %>% 
  dplyr::mutate(HeightMax = as.numeric(HeightMax)) %>% 
  dplyr::mutate(DiameterMax = as.numeric(DiameterMax)) 

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

Filter mounds

mnds <- features %>% 
  dplyr::filter(Type == "Burial Mound" | Type == "Extinct Burial Mound" | Type == "Uncertain Mound")
mnds

## Simple feature collection with 1304 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 420515 ymin: 4647591 xmax: 499612.3 ymax: 4722135
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1,304 × 19
##     TRAP Source             Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##  * <dbl> <chr>              <chr> <chr>   <chr>    <dbl>   <dbl> <fct>   <chr>  
##  1  9305 <NA>               Buri… Scrub   Scrub       25     1.8 3       <NA>   
##  2  9411 Survey             Buri… Scrub   Scrub       15     1.5 2       Nodata 
##  3  9416 Legacy verificati… Buri… Annual… Annua…      20     2   3       Agricu…
##  4  9417 Legacy verificati… Buri… Pasture Pastu…      15     2   4       Looting
##  5  9409 Legacy verificati… Buri… Annual… Pastu…      30     2   2       Agricu…
##  6  9314 Legacy verificati… Buri… Pasture Pastu…      30     2   5       Looting
##  7  9315 Legacy verificati… Buri… Pasture Pastu…      20     1   2       Looting
##  8  9415 Legacy verificati… Buri… Annual… Annua…      25     1   5       Agricu…
##  9  9418 Legacy verificati… Buri… Pasture Pastu…      50     6   2       Looting
## 10  9319 Legacy verificati… Buri… Annual… Pastu…      15     1   3       Agricu…
## # … with 1,294 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## #   distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## #   prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## #   variable names ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition,
## #   ⁵​PrincipalSourceOfImpact

Create buffers of 200 m

mnds_buff200 <- st_buffer(mnds,200)

Explore clusters

How many mounds have one or more neighbor(s) within 100 m? Let’s calculate the number of neighbors via a sdgp matrix with st_intersects() and then sum() and filter() help us determine their identity. We can review their location with mapview()

library(mapview)
mapview(mnds_buff200) + mapview(mnds)

test <- st_intersects(mnds_buff200, mnds_buff200$geometry)
test0 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>0 # intersects at least one - everything as each intersects itself
test1 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>1 # 848
test2 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>3 # 361 mounds in in clusters of 2+ mounds
test5 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>6 # 147 features 95-97
test6 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>7 # 119
test7 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>8 # 87
test10 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>11 # 30
test11 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>12 # none
test21 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>20 # none

# How many clusters of 2 or more mounds are there? 
sum(test2)

## [1] 361

# How many clusters of 10?
sum(test7)

## [1] 87

# Which mounds are in clusters of 10?
mnds %>% filter(test10)

## Simple feature collection with 30 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 467542.2 ymin: 4683521 xmax: 478932.4 ymax: 4721394
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 30 × 19
##     TRAP Source             Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##  * <dbl> <chr>              <chr> <chr>   <chr>    <dbl>   <dbl> <fct>   <chr>  
##  1  8022 RS:SITE            Buri… Annual… Annua…      20     0.7 5       Agricu…
##  2  8024 RS:SITE            Exti… Annual… Annua…      20     0.5 5       Agricu…
##  3  8023 RS:SITE            Exti… Annual… Annua…      25     0.1 4       Agricu…
##  4  8025 RS:SITE            Exti… Annual… Annua…      20     0.5 5       Agricu…
##  5  8026 RS:FNEG            Exti… Annual… Annua…      20     0.3 5       Agricu…
##  6  8027 RS:FNEG            Buri… Annual… Annua…      30     1   4       Agricu…
##  7  8028 RS:FNEG            Buri… Annual… Annua…      40     1   4       Agricu…
##  8  8021 RS:FNEG            Exti… Annual… Annua…      20     0.3 5       Agricu…
##  9  8360 Survey             Buri… Annual… Annua…      40     1.4 3       Agricu…
## 10  8365 Legacy verificati… Buri… Annual… Annua…      35     1.1 3       Agricu…
## # … with 20 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## #   distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## #   prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## #   variable names ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition,
## #   ⁵​PrincipalSourceOfImpact

Visualise clusters among the Yambol mounds

# View the clusters of 5 to 10: 
# We add the residual mounds so that we see all the neighbors as not all get tagged!
mnds %>% filter(test5) %>% 
  mapview() + mapview(mnds, cex = 0.1, zcol = "Type")

At least 10 mounds come out as having 10 neighbors within 100m, however, each of these clusters contains duplicates - mounds which were registered repeatedly under different numbers in successive seasons. Not all of the mounds inside these clusters are marked as members of 10-mound-cluster, as only the central ones have all 10 within their buffer. The mounds at the far edges of the cluster do not get flagged as their buffers do not contain all the cluster fellows.

Furthermore, at least two clusters contain duplicates: - In the cluster W of Karavelovo, the three duplicates are 8022, 8023, 8024, and 8028 which were rudimentarily recorded in 2009 as part of remote sensing ground truthing, and again in 2010 as 9594, 9593, and 9595 and 9592 (from 2010+). IN this case 2009 can be eliminated. This means this cluster has 4 mounds less than counted. - There is another cluster near the quarry west of the Mogila village, with the following five duplicates: 8346, 8347, 8348, 8349, 8350 from 20XX to 9222, 9223, 9224, 9225, 9226. These need to be reviewed.

Check duplicates

# Sanity check on the duplicates in the other clusters
mnds[95,]

## Simple feature collection with 1 feature and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 464359.3 ymin: 4706914 xmax: 464359.3 ymax: 4706914
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1 × 19
##    TRAP Source              Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##   <dbl> <chr>               <chr> <chr>   <chr>    <dbl>   <dbl> <fct>   <chr>  
## 1  9222 Legacy verification Buri… Pasture Pastu…      20       2 3       Looting
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## #   elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## #   TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## #   ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition, ⁵​PrincipalSourceOfImpact

test[[95]] # buffer around point 95 has allegedly over 6 neighbours

##  [1]  95  96  97 347 348 609 610 613 614 616 617

#[1]95  96  97 347 348 610 613 614 616 617
mapview(mnds[c(95:97,347:348,613:614,616:617),])

# Duplicates - Mogila cluster!
dplyr::distinct(mnds)

## Simple feature collection with 1304 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 420515 ymin: 4647591 xmax: 499612.3 ymax: 4722135
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1,304 × 19
##     TRAP Source             Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##    <dbl> <chr>              <chr> <chr>   <chr>    <dbl>   <dbl> <fct>   <chr>  
##  1  9305 <NA>               Buri… Scrub   Scrub       25     1.8 3       <NA>   
##  2  9411 Survey             Buri… Scrub   Scrub       15     1.5 2       Nodata 
##  3  9416 Legacy verificati… Buri… Annual… Annua…      20     2   3       Agricu…
##  4  9417 Legacy verificati… Buri… Pasture Pastu…      15     2   4       Looting
##  5  9409 Legacy verificati… Buri… Annual… Pastu…      30     2   2       Agricu…
##  6  9314 Legacy verificati… Buri… Pasture Pastu…      30     2   5       Looting
##  7  9315 Legacy verificati… Buri… Pasture Pastu…      20     1   2       Looting
##  8  9415 Legacy verificati… Buri… Annual… Annua…      25     1   5       Agricu…
##  9  9418 Legacy verificati… Buri… Pasture Pastu…      50     6   2       Looting
## 10  9319 Legacy verificati… Buri… Annual… Pastu…      15     1   3       Agricu…
## # … with 1,294 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## #   distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## #   prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## #   variable names ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition,
## #   ⁵​PrincipalSourceOfImpact

duplicates <- lengths(st_equals(mnds))>1 # produces zero duplicates, as they do not have identical locations, but are withing 10-20m of one another. A spatial filter on close proximity is needed

# Close mounds
mnds[c(95:97,347:348,613:614,616:617),]

## Simple feature collection with 9 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 464228.7 ymin: 4706914 xmax: 464404.7 ymax: 4706973
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 9 × 19
##    TRAP Source              Type  LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
##   <dbl> <chr>               <chr> <chr>   <chr>    <dbl>   <dbl> <fct>   <chr>  
## 1  9222 Legacy verification Buri… Pasture Pastu…      20    2    3       Looting
## 2  9225 Legacy verification Buri… Pasture Pastu…      15    1    2       None   
## 3  9226 Legacy verification Buri… Pasture Pastu…      20    1.5  3       Looting
## 4  9223 Survey              Buri… Pasture Pastu…      15    1.5  4       Looting
## 5  9224 Survey              Buri… Pasture Pastu…      20    1.5  3       None   
## 6  8347 Legacy verification Buri… Pasture Pastu…      17    1.2  1       Post-d…
## 7  8348 Legacy verification Buri… Pasture Pastu…      15    1.5  2       Looting
## 8  8349 Legacy verification Buri… Pasture Pastu…      17    1.75 2       Looting
## 9  8350 Legacy verification Buri… Pasture Pastu…      18    1.2  2       Post-d…
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## #   elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## #   TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## #   ¹​LU_Around, ²​DiameterMax, ³​HeightMax, ⁴​Condition, ⁵​PrincipalSourceOfImpact

Voronoi polygons

There is ast_voronoi() function but when using it on points, one must st_combine() the points first (as with convex hull) so as to create a mesh covering all the points.

# Create Voronyi geometries
vor <- st_voronoi(st_combine(mnds))

# st_voronoi returns a GEOMETRYCOLLECTION,
#  some plotting methods can't use a GEOMETRYCOLLECTION.
#  this returns polygons instead
mnd_vor_poly <- st_collection_extract(vor)

# Crop the polygons to Yambol region boundary
vor <- st_intersection(Y_region, mnd_vor_poly)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# Visualize
mapview(vor)+ mapview(mnds %>% filter(test5)) + mapview(mnds, cex = 0.1)

Nice overview of where some of the medium sized clusters are. Next step: mark the clusters by 1-10 depending on how many neighbors each mound has and produce a colour-coded version so one can see the clusters better.